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1  . 

Divide  and  conquer  (DC)  methods  have  been  developed  for  the  symmetric 
e i genprob 1 em ,  see  Cuppen  [Cu] ,  Dongarra  and  Sorenson  [DS] ,  and 
Krishnakumar  and  Morf  [KM] ,  and  have  for  these  problems  been  shown  to  be 
efficient  on  parallel  computers  and  competitive  on  single  processor 
machines  [DS]  ,  [Cu]  ,  [KM]  .  The  DC  method  has  also  been  appl ied  successful ly 
to  the  computation  of  singular  values  by  Jessup  and  Sorensen  [JS] .  In  the 
present  paper  we  describe  a  DC  method  for  the  unitary  e igenprobl em ,  and  we 
also  discuss  the  simplifications  that  arise  for  real  orthogonal  matrices. 

Let  H  €  Cnxn  be  unitary.  Then  H  is  unitari ly  similar  to  a  upper 
Hessenberg  matrix  with  real-valued  non-negative  subdiagonal  elements.  If 
a  subdiagonal  element  vanishes,  then  the  eigonproblem  splits  into 
e  i  genprobl  err.s  for  smaller  upper  Hessenberg  matrices.  We  therefore  may 
assume  that  H  is  a  upper  Hessenberg  matrix  with  pos i t i ve  subdiagonal 
elements.  Then  all  eigenvalues  of  H  are  simple.  It  is  easily  seen  that  H 
can  be  written  as  a  product  of  n  Givens  reflectors  Gj  G  Cnxn, 

H  =  H(7j,72»-..i')n):=Gx^2,,'^n-l^'n’  (  1  •  1  ) 

where,  for  1  <  k  <  n, 


->k 

17  k  f  k 


7k  G  C,  trk  G  R,  <7k  >  0, 
I  Tk  1 2  +  ‘’’k2  =  1  > 


and 


Gn  :  = 


1  n-1 


“  7  n 


!  7n  G  C,  |  7n  I  —  1 


(1.2a) 


(1.2b) 


2 


Here  Ij  denotes  the  jxj  identity  matrix.  The  ,  1  <  j  <  n,  are  the  so- 

called  Schu r  paramete rs  of  H,  and  7j  denotes  the  complex  conjugate  of  ij . 
The  tTj ,  1  <  j  <  n,  are  said  to  be  com  pi  ementary  paramete  rs  of  H,  and  are 
the  subdiagonal  elements  of  H. 

The  DC  method  described  uses  the  product  representation  (1.1)  of  H, 
the  so-called  Schu r  paramet r i c  f o rm  of  H.  An  application  of  particular 
interest  to  us  is  the  computation  of  Pisarenko  frequency  estimates  for  a 
random  stationary  stochastic  process,  see  below.  In  this  application  H 
defined  by  its  Schur  parametric  form.  The  determination  of  Gaussian 
quadrature  rules  on  the  unit  circle,  so-called  Gauss-Szego  quarature 
rules,  also  gives  rise  to  unitary  (or  real  orthogonal)  matrices  in  Schur 
parametric  form,  see  Section  5.  When  the  Schur  parameters  are  not 
explicitly  known,  they  can  be  computed  from 


i  s 


ll 


-II 


ll 


Ti 


°2  C1  H>ij.  J  =  2’3 . . 


where  G^  denotes  the  conjugate  transpose  of  Gk  ,  and  Mjj  denotes  the  element 
(j,j)  of  a  matrix  M  G  CnXn, 

The  DC  method  is  most  easily  described  for  H  G  RnXn  orthogonal.  Then 


Hi 

U 

H  =  Gt  G2  . 

■  •  Gs  Gs+1  . 

.  .  Gn  = : 

I  n-s 

Gs 

h2 

where  Hj  G  Rsxs  and  H2  G  s)x(n  s)  are  orthogonal.  The  Givens  reflector 
Gs  G  Rnxn,  s  <  n,  can  be  written  as  a  Householder  transformation 

Gs  =  I  -  2wwh  (1.4) 


where  : 


3 


w  :  = 

esws  +  es^_1ws^_1  €  Rn  , 

(1.5) 

ws  :  = 

2  1/2 ( 1  +  ^s)l/2  , 

(1 .6a) 

“s+i 

:=  -2  1/2  ( 1  -  ^)l/2  . 

(1.6b) 

Throughout  this  paper  ej  denotes  the  jth  column  of  an  identity  matrix  of 
appropriate  order.  By  (1.3)-(1.4),  H  is  orthogonally  similar  to 


H' 


(I  -  2wwH ) 


H  -  2HwwH  . 


(1.7) 


This  is  one  step  of  the  DC  method  for  orthogonal  matrices:  The  eigen¬ 

values,  and  if  so  desired  the  eigenvectors,  of  Hj  and  H2  are  computed 
first.  Hr  is  a  rank-one  modification  of  H,  and  the  eigenvalues  of  H  are 
computed  as  the  zeros  on  the  unit  circle  of  a  rational  function,  whose 
poles  are  the  eigenvalues  of  H* . 

Section  2  describes  the  DC  method  for  unitary  matrices.  In  Section  3 
we  show  some  results  on  the  orthogonality  of  the  eigenvectors  and  on  the 
location  of  the  eigenvalues.  These  results  are  analogous  to  bounds 
presented  by  Dongarra  and  Sorensen  [DS]  for  the  DC  method  for  symmetric 
matrices.  Section  4  discusses  simplifications  that  can  be  made  when 
H  is  real  and  orthogonal,  and  also  considers  some  computational  details. 
Computed  examples  for  the  orthogonal  eigenvalue  problem  are  presented  in 
Sect  ion  5 . 

An  outline  of  a  unitary  DC  method  with  convergence  results  for  the 
root-finder  has  previously  been  presented  in  [GR] .  The  splitting  into 
subproblems  is  done  differently  in  the  present  paper.  A  related  DC  method 
is  described  by  Arbenz  and  Golub  [AGJ .  Cybenko  [CyJ  reduces  the 
orthogonal  eigenproblem  to  an  eigenproblem  for  a  symmetric  tridiagonal 


4 


matrix.  The  orthogonal  eigenproblem  is  in  [AGR1]  solved  by  solving 
singular  value  problems  for  certain  bidiagonal  matrices,  and  a  QR 
algorithm  for  unitary  matrices  is  presented  in  [Grl] .  A  comparison  with 
respect  to  accuracy  and  speed  of  these  methods  still  remains  to  be  done. 
Here  we  only  note  that  DC  methods  are  suitable  for  implementation  on 
parallel  computers,  see  [DS] , [KM] ,  and  Section  2.  Some  of  the  schemes 
mentioned  transform  the  orthogonal  eigenvalue  problem  to  a  symmetric  one. 
Tiie  real  eigenvalues  of  the  latter  problem  are  then  mapped  to  the  unit 
circle  to  yield  the  eigenvalues  of  the  orthogonal  eigenproblem.  Tiie 
mapping  from  the  interval  to  the  unit  circle  may  be  sensitive  to 
perturbations  of  arguments  near  the  end  points  of  the  interval ,  and  it  may 
therefore  be  difficult  to  determine  eigenvalues  close  to  ±1  accurately 
with  these  schemes. 

Pisarenko  [Pi]  proposed  a  method  for  decomposing  a  random  stationary- 
stochastic  process  {xm}^_g,  xm  €  R,  into  a  sum  of  harmonics  in  white 
no i se ,  i . e . , 

P 

Xm  =  E  Qp  cos  +  0  A  +  ym  ,  m  >  0 ,  (1.8) 

e=i  it 

where  the  6 ^  are  arbitrary'  phase  shifts  and  {ym}^’_g  is  a  zero  mean  white 
noise  process  with  variance  <x2 .  The  4>  g  are  called  P  i  sarenko  frequency' 
est i mates .  Assume  for  simplicity  that  p  is  the  number  of  distinct 
harmonics  in  the  ‘signal’  {xm}^=()  is  known,  and  that  0  <  <j>^  <  tt  for 
1  <  ^  <  P-  Then  the  < i>^  can  be  determined  as  follows.  Form  the 

(2p+ 1 ) x (2p+ 1 )  Toeplitz  covariance  matrix  M  for  the  signal  {xm)^_jj,  and 
compute  its  least  eigenvalue  Amjn  .  Then  Amjn  =  a2,  see  [Pi].  Let  be 

the  Schur  parameters  associated  with  the  Toeplitz  matrix  M-A  (  1 .  They  can 
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be  determined  from  the  Szego  recursions  (Levinson’s  algorithm),  see  e.g. 

[AGR2]  .  From  our  assumptions  it  follows  that  M-AminI  is  singular,  but 

leading  principal  submatrices  not  identical  with  M-AminI  are  not. 

Therefore,  - 1  <  7j  <  1  for  1  <  j  <  2p,  and  y2p  £  {-1,1}.  By  (1.1)-(1.2)  it 

2P 

follows  that  the  Schur  parameters  {'h}T.  define  an  orthogonal  matrix  H  € 
R2pX2p  with  distinct  eigenvalues  Enumerate  the  eigenvalues  so  that 

those  with  Im(Aj)  >  0  have  smaller  index  than  the  eigenvalues  with  Im(Aj)  < 
U.  Then  the  Pisarenko  frequency  estimates  are  given  by 

:=  a  r  g  ( A  j )  ,  1  <  j  <  p  . 

The  coefficients  Oj  of  (1.8)  are  two  times  the  weights  belonging  to  the 
Gauss-Szego  quadrature  rule  with  abscissas  Aj  ,  1  <  j  <  p.  For  details  see 

[AGR2] ,  where  also  references  to  related  work  can  be  found.  The  unitary 
DC  method  yields  the  Gauss-Szego  weights  with  no  extra  computational 
effort  when  computing  the  eigenvalues  Aj  .  Gauss-Szego  quadrature  is 
discussed  in  Example  5.1  of  Section  5. 
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The  unitary  e  i  genprob  1  etn 


2  . 

In  this  section  we  describe  a  divide  and  conquer  method  for  unitary 
right  Hessenberg  matrices  with  positive  subdiagonal  elements.  First  we 
need  to  generalize  the  splitting  (1.3)-(1.7)  to  Givens  reflectors  with 
complex-valued  Schur  parameters.  This  is  accomplished  by  noting  that  the 
Gs  defined  by  (1.2a)  are  diagonally  unitarily  equivalent  with  a  real 
Householder  transformation.  Introduce 


Is  :  — 


Is/  I  7s  I  7s  ^  0 

1  ,  Is  =  0  • 


Then  |  I  =  1'  and  for  Gs  defined  by  (1.2a)  we  obtain 


- 1 

+ 

C/l  — . 

_ I 

Cs 

Is 

7s 

In-S 

I  n-s-1 

IS-1 

_  1 7  s  1 

*s 

^s 

1  7s  1 

- 

^  n-s-1 

,  similarly 

to 

(1 

•  5) - ( 1 .6) 

'  ’  =  esws 

(2.1) 


ws  :=  2'1/2(1  +  |  7s  I  )  1/2 


Js+l 


:=  ~  2  1/2  (  1  -  |7s|) 


l/2 


(2.2a) 

(2.2b) 


Substitution  of  (2.1) 


II 


I 


n-s 


(i 


into  (1.1) 

h  r ** 

2wwH) 


yields 


ll2 


(2.3) 


with 


Hi  H  ( 7 1 , 7  2  >  •  •  •  i  7s-i  >  ~~ls  )  > 

H2  :=  H  (7sSs+i  »  TsSs+2  »  •  •  •  »  TsSn)  • 


A  unitary  similarity  transform  of  (2.3)  yields,  analogously  with  (1.7), 

H'  :  = 

Let 

Hk  =  Vk  Ak  W»  ,  k  =  1,  2,  (2.5) 


(I  -  2wwH)  =:  H  -  2HwwH 


(2.4) 


be  spectral  resolutions,  i  .e.  ,  the  Wk  are  unitary  and  the  Ak  are  diagonal  . 
Then  11  has  the  spectral  resolution  H  =  W  A  WH ,  where 


W 


6  i  aS  ( ^1  * ^2  ’  •  ’  '  >  ) 


(2.6) 


with  |  Ak  |  =  1  for  1  <  k  <  n. 

We  are  in  a  position  to  describe  how  the  spectrum  of  H  can  be  obtained 
from  A,  the  last  row  of  W!  and  the  first  row  of  W2 •  Introduce  the 
characteristic  polynomial 


X  (A)  :=  det(AI-H)  =  det(AI-Il')  =  det  ( AI -H+2Hwwh) 

=  det  (AI  -H)  det(  I +2  (  A I  -ll  )'1FlwwH) 

=  V»(  A)(l  +  2wH(AI-fi)‘1fiw) 

=  1>{  A)  (  1  +  2whW  (  A  I  -  A  ) A  WHw)  , 


where  H;  is  defined  by  (2.4)  ,  W  and  A  by  (2.6)  ,  w  by  (2.2)  and  V'(A)  :  = 

det(A[-H]).  Let  z  =  [CDjLj  be  given  by 


z 


WHw 


V^esws 
W2  elws+l 


(2.7) 


and  define  the  spectral  function 


*(A) 


X^A)  _ 

V-'(A)  ~ 


+  2zH  (A  I  -  A  )_1  Az 


n 

1  +  2  £ 
j=l 


=  i  icji 

j=i 


2 


A  -f  Aj 
A  -  Aj 


(2.8) 


H 

where  we  have  used  that  z  z  =  1 .  Let 

8j  :=  arg(Aj)  ,  8  :=  arg(A)  ,  0  <  6j ,  6  <  2ir 


Then,  with  i  :=  4-1 , 


n  /9-  - 

<HA)  =  i  !  Cj  1 2  cot^-b-^ — j  =:  i<l>(0)  , 


Q'C#)  =  5  I  Cj !  2/s  i  n2(^-2 — )  >  5  zHz  =  2 


(2.9) 


(2.10) 


We  may  assume  that  the  are  distinct  and  that  all  <t  ^  0,  because 
otherwise  we  can  make  these  conditions  hold  by  deflation,  see  below.  Let 
6-*  &  [0,2*[,  1  <  j  <  n,  denote  the  zeros  of  $(0).  Then  the  eigenvalues  of 

H*  and  of  H  are  given  by  Xj  :  =  exp(i0jf),  1  <  j  <  n.  The  sets  and 

{$/  }j=  ]  strictly  interlace. 

We  describe  a  rootfinder  for  <J>(0).  By  the  inequality  (2.10),  the 
zeros  of  <J>(0)  can  be  determined  accurately.  We  may  assume  that 
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0  <  6 1  <  $2  <  •  •  •  <  #rt  <  2ir  and  that  0^°\  our  initial  approximation  of  a 


zero  of  $(0)  ,  satisfies  0n-2ir  <  6 


(0) 


By  the  strict  interlacing  of  the 


sets  { 0 j } J1—  |  and  ,  $(0)  has  precisely  one  zero,  denoted  6 x  .  in  the 

open  interval  ]0n-2jr,  9X  [ .  Assume  for  the  moment  that 


$(0(o))  <  0  , 


(2.11) 


and  introduce 


<£(#)  :=  p  +  a  cot 


r-h r -)  ■ 


(2.12) 


The  coefficients  p  and  a  are  determined  by  osculatory  interpolation,  i.e.. 


4>(0(o))  =  4>(0(o))  .  4>'(0(o))  =  <f>'(0(o))  , 


(2.13) 


which  vie  Ids 


(p  =  <1 >(0(o))  -  <D'(0(o))  sin(Ji  -  0(o))  , 

Iff  =  2<t>'(0(o))  s  i  n2^  —  ~r26 - . 

Ilie  zero  0  ^  of  ^(O)  in  ']0n-2v,6x[  is  our  next  approximation  of  0X  .  New 

approximations  0(j+1*  of  $x  are  computed  from  0(j) ,  j  >  1,  in  a  similar 
fashion.  Ilie  sequence  satisfies  0^  <  9X  for  j  >  0,  and  converges 

monoton i cal  1 y  and  quad  rat i ca 1 1 y  to  0 x  as  j  increases,  see  [GR]  for  a 
proof . 

If  instead  of  (2.11)  we  have 


$(0(O))  >  0  , 


then  we  replace  (2.12)  by 


4>(0)  =  p  +  a  cot(^_- — 


(2.14) 


(2.15) 
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of  <t>(0)  in  the  open 


in  (2.13).  This  defines  p  and  a .  The  zero 
interval  ]Sn-2i,^1[  is  our  next  approximation  of  .  New  approximations 
^(j+i)  0f  are  computed  front  0(j)  for  j  >  1  in  a  similar  fashion.  The 

sequence  {0^  Ij^o  satisfies  0^  >  6 ^  for  j  >0,  and  converges  monoton  i  cal  1  y 
and  quad  rat  i  ca  1  1  y  to  6-^ ,  see  [GR]  .  In  the  implementation  used  to  generate 
the  computed  examples  of  Section  5,  the  iterations  are  carried  out  until 
$(0^J"*"1^)  >  4>(0^  <  <$(0^  1^)  .  The  value  0^  is  accepted  as  an  approximate 

root  of  $  (0)  =  0  . 


F  rom 


r  10  ' 

A  :  =  d  i  ag  [e 


i0. 


(2.16) 


and  the  spectral  resolutions  (2.5)  of  1^  and  H2 ,  we  can  now  compute  the 
spectral  resolution  of  H: 


II  =  WAWh  ,  VHV  =  I  . 


(2.17) 


By  (2. 3) -(2. 4),  we  can  for  some  vector  u  G  Cn  express  H  as 


11  =  II  -  2uu 


~  H 


u  :  = 


l^eju-s 


elu's+l 


G  Cn  . 


Let  A  exp(i0r)  and  v  G  Cn  form  an  eigenpair  of  H,  i.e.  Ilv  =  vA .  Then 


(li  -  2 u u H )  v  =  vA  , 


or,  equivalently, 


(II  -  I A )  v  =  Do  ,  q  :=  2uHv 


Ihis  shows  that  any  normalized  eigenvector  v  of  II  associated  with  the 
eigenvalue  A  is  a  normal ization  of 


where  Wk  and  Ak  are  given  by  (2.5).  Let  ||  ||2  denote  the  Euclidean  vector 
and  matrix  norms.  From  ||Wj||2  =  ||W2||2  =  |J  A  j  ||2  =  1  and  (2.6),  (2.7),  (2.10), 

(2.18)  it  Follows  that 


«(A)  :=  l|v'j|2  =  ||  (A- 1  A)'1  x2  j|2 


n 

=  (£  nr 

J  =  1  I  Aj 


y'2  =  (i*'(»))‘/2  >  i 


(2.19) 


We  choose 


v 


A 


W1(I-Aj,A)-1wJlesws 

W’2(A2-IA)  'w^ejU^j 


/<(A)  , 


(2.20) 


and  note  that  the  lower  bound  (2.19)  for  6(A)  indicates  that  severe 
cancellation  of  significant  digits  does  not  take  place  in  the  computation 
of  vA  by  (2.20) . 

By  (2.7)  ,  we  only  require  the  last  row  of  Wx  and  the  first  row  of  W'2 
(as  well  as  A)  in  order  to  determine  the  spectral  function  (2.8).  Hence, 
if  we  do  not  desire  the  eigenvectors,  then  it  suffices  to  determine  the 
first  and  last  rows  of  W /  in  order  to  be  able  to  compute  the  spectrum  of 
larger  problems.  We  call  the  triplet  {A,e|*W,en  W}  the  part  i  al  spectral 
reso 1 u t ion  of  H.  The  first  and  last  elements  of  v^  can  easily  be 
determined  from 


felIvA  =  e,1IW1(I-AfA)1wJ,esu;s/6(A)  , 
le”vA  =  e»sW2(A2-IA)-1w“e1u;s+1/6(A) 
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We  may  assume  that  the  columns  of  W  are  such  that  all  components  of  the 
vector  WHex  are  real  and  positive.  Then  e^WHej  is  the  weight  correspond¬ 
ing  to  the  node  exp(i0k/)  in  the  Gauss-Szego  quadrature  rule  with  nodes 
{exp(i^j,)}j'_j,  see  [Gr2]  and  Example  5.1. 

We  assumed  above  all  components  <(g  of  z  to  be  non-vanishing  and  all 
eigenvalues  Aj  of  H  to  be  distinct.  These  conditions  can  be  made  to  hold 
true  by  def 1 at  ion.  Our  discussion  follows  Dongarra  and  Sorensen  [DS]  . 
First  assume  that  (g  vanishes.  By  (2.4)-(2.7), 

WhH'V  =  A  -  2WhHwwhW  =  A  (I  -  2zzH)  ,  (2.22) 

1  i 

and  from  eg  z  =  0  it  follows  that 

A  (  I  -  2zzH)eg  =  Ae(  =  Ageg,  Ag  :=  e{HAe£  .  (2.23) 

Substituting  (2.23)  into  (2.22)  yields 
HfWeg  =  WegAg  , 
and  therefore 


Ihus  if  Q  =  0,  then  we  can  determine  an  eigenpair  of  H  without  explicitly 
computing  a  zero  of  $(0)  and  without  using  (2.20). 

For  a  general  z  €  Cn ,  with  zHz  =  1 ,  we  obtain 

||A(I-2zzH)eg||2  =  2  |  C£  I  , 


and  we  accept  {Ag,eg}  as  an  eigenpair  of  A  if 


2  |  C£l  <  <1 


(2.24) 


for  some  small  constant.  ^ . 

Assume  that  z  =  [CjjjLj  with  2  |  <7  |  >  for  al  1  j  .  We  may  now  be  able  to 

deflate  due  to  close  eigenvalues.  Let  A  =  d  iag  [Ax  ,  A2  ,  .  .  .  ,  A„]  with 
Ax  «  A2 ,  and  choose  the  Givens  reflector 


G 


-7  <r 
a  7 


I 


n-2 


6 


Cnxn  , 


(2.25) 


so  that 


Gz  =  C(  I  Cl  1  2  +  K2I2)1/2>0’<3»C4.  •  •  •  ><n]T  , 


i  .  e  . 

( a  s=  K2|/(|fil2  +  K2I2)1/2  , 

U  :=  ~<i  |-^-|/(  I  Cl  1 2  +  K2|2)l/2  • 

We  accept  {h<r2  +  72l7|2>GHe2}  as  an  (approximate)  eigenpair  of 
A ( I  -  2zzH)  if 

|7 »(>i  -  >2)  I  <  f2,  (2.26) 

for  some  smal 1  constant  e2 ,  because 

||A(I-2zzH)GHe2  -  (Ax(72  +  72  I  7 1  2)GHe2||2  =  |  7^  (7a  -  72)  I  • 

If  7j  =  72  then  we  have  determined  an  eigenpair  exactly.  In  case 

Aj  ^  A2  we  note  that  |7«r(Aj  -  A2)  |  <  g  |  Aj  -  A2  |  ,  and,  moreover,  if  1 7 1  «  0 

or  1 7 1  ss  1  then  ^^(Aj  -  A2)  |  <<  |  Ax  -  A2  |  .  Hence,  inequality  (2.26)  may 

be  satisfied,  even  if  | Ax  -  A2  |  >  e2 .  Assume  that  (2.26)  is  valid.  Then  A 

is  replaced  by 
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A  :=  d  i  ag  [Aj  ("y2)  +  A2<t2  ,  A3  ,  A4  ,  .  .  .  ,  An]  €  C(n-lMn*1' 

and  if  A  has  close  eigenvalues,  then  deflation  is  repeated. 

The  unitary  DC  method  can  be  used  in  two  ways.  One  approach  is  to 
divide  the  original  e igenproblem ,  as  well  as  subproblems  so  obtained, 
until  only  trivial  e i genprob 1 ems  of  orders  two  and  one  remain.  These 
small  e igenprob 1 ems  are  solved  analytically.  From  the  solutions  of  small 
e igenproblems ,  the  solutions  of  e igenprobl ems  of  larger  size  are  computed, 
and  this  step  is  repeated  until  the  solution  of  the  original  eigenproblem 
has  been  determined.  This  approach  is  used  in  the  numerical  examples  of 
Section  5. 

An  alternative  approach  is  to  use  the  DC  technique  to  generate  just  a 
few  sube i genprob 1 ems ,  each  of  which  can  be  solved  independently  by  some 
other  numerical  scheme,  such  as  the  unitary  GR  method  [Grl] ,  or  the  scheme 
in  [AGR1] ,  in  case  the  matrix  is  real  orthogonal . 

We  conclude  this  section  with  some  bounds  of  the  computational 
complexity  of  the  unitary  DC  method.  Assume  that  H  G  Cnxn  is  given  in 
Schur  parametric  form  (1.1)  with  positive  subdiagonal  elements  (Tj .  Let  n  r= 
2  for  some  positive  integer  £,  and  subdivide  the  given  eigenproblem  until 
2  e i genprobl ems  for  2x2  matrics  are  obtained.  The  latter  e i genprob 1 ems 
are  solved  analytically.  We  assume  that  the  number  of  iterations  required 
by  the  rootfinder  for  4>(0)  can  be  bounded  independently  of  n. 

Let  first  n  independent  processors  be  available.  The  reduction  of  the 
original  eigenproblem  for  H  to  ^  e igenproblems  for  2x2  matrices  can  be 
carried  out  in  t1  :=  O(log22)  time  steps.  This  computation  only  requires 
the  determination  of  the  Schur  parameters  for  the  unitary  matrices  of  the 
smaller  e i genprob 1 ems ,  see  (2.3).  Let  the  Schur  parameters  for  all  ^ 
unitary  2x2  matrices  be  known.  The  spectral  resolution  of  all  these 
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Assume  "that  the 


matrices  can  be  computed  in  t2  :=  0(1)  time  steps, 
partial  spectral  resolutions  (2.21)  of  all  2^*"*"1  unitary  2*  1  x  2*  1  matrices 
are  known  for  some  j  €  [2 , £]  .  In  order  to  compute  the  partial  spectral 

resolution  of  all  2  unitary  2s  x  21  matrices,  we  have  to  compute  2*  zeros 

C-i 

of  each  of  the  2  functions  $(0)  ,  see  (2.8).  Hence,  a  total  number  of  n 
zeros  have  to  be  computed,  and  we  use  one  processor  to  determine  each  one. 
Each  function  $(0)  has  2*  terms,  and  can  therefore  be  evaluated  in  0(2J) 
time  steps  for  each  value  of  6.  Hence,  we  can  determine  all  eigenvalues 
of  all  2^*  unitary  2*  x  2*  matrices  in  tg^  :=  0(2*)  t  ime  steps.  For  each 
eigenvalue  we  compute  the  first  and  last  elements  of  the  corresponding 
eigenvector  from  (2.21).  The  first  and  last  element  of  one  eigenvector 
can  be  determined  by  one  processor  in  0(2*)  time  steps.  These 
computations  have  to  be  carried  out  for  n  eigenvectors  by  n  processors  and 
therefore  require  t^  =  0(2*)  time  steps.  Hence,  the  number  of  time  steps 
required  to  determine  the  partial  spectral  resolution  of  H  by  n  processors 
i  s 

l  (j)  ^  (i) 

tx  +  t2  +  £  ti  +  £  t  =  0  (n )  .  (2.27) 

j=2  j=2 

Now  let  n2  independent  processors  be  available,  and  assume  that  the 
partial  spectral  resolutions  of  all  2^*+1  unitary  2*"1  x  2*"1  matrices  are 
known  for  some  j  £  [2 , £]  .  We  have  now  n  processors  available  for  each 

evaluation  of  each  of  the  2  functions  $((?).  Each  of  these  functions  $(0) 
has  21  terms  and  can  for  each  value  of  9  be  evaluated  in  0(log22*)  time 
steps.  Hence,  we  can  compute  al  1  eigenvalues  of  all  2^  *  uni  tary  2*  x  2* 
matrices  in  tg  ^  :=  0(log22*)  time  steps.  The  first  and  last  elements  of 
each  eigenvector  of  each  of  the  2^"*  unitary  2*  x  2*  matrices  can  be 


determined  in  t ^  :=  0  ( 1  og2  2? )  time  steps,  by  using  n  processors  to 

compute  each  sum  with  2*  <  n  terms.  The  initial  determination  of  the  Q 
unitary  2x2  matrices  and  their  spectral  resolutions  cannot  be  sped  up 
essentially  by  using  more  than  n  processors,  and  requires  0(tj)  +  0(t2) 
time  steps.  Hence,  the  number  of  time  steps  required  in  order  to  compute 
the  partial  spectral  resolution  of  H  by  n2  processors  is 


0(t!) 


+  0(t2) 


t 

+  £ 

1=2 


0(1  og2  n )  +  £  0(j) 
1=2 


0(1  og^  n  ) 


(2.28) 


The  time  complexities  (2 . 27) - (2 . 28)  suggest  that  the  unitary  DC  method 
presented  could  be  attractive  for  use  in  real-time  signal  processing 
app 1 i cat  ions. 
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3 .  Some  properties  of  the  unitary  DC  method 


We  show  some  properties  of  the  eigenvectors  of  H  and  zeros  of  Q>(6). 
Analogous  results  have  previously  been  obtained  by  Dongarra  and  Sorensen 
[DS ,  Lemmas  4.2,  4.6  and  4.7]  for  the  DC  method  for  the  eigenproblem  for 
symmetric  tridiagonal  matrices.  The  following  formulas  are  used  in 
several  of  the  proofs: 


H  A)  =  1+2  t  I  Cj  I  2 

j  =  l 

*'(A)  =  -2  t  ICjl2  ~ 

j  =  l  (A 

*'(8)  =  <t>\  A)A  . 


(3.1) 

(3.2) 

(3.3) 


Lemma  3 . 1  Let  Aj/i  G  C,  j  A  |  =  |/i|  =1  and  A  ^  Assume  that  g  { Aj }  j  , 

where  A  =  d  i  ag  [Aj , A2 ,  .  .  .  , An]  .  Let  v^  and  v^  be  defined  by  (2.20).  Then 


I  v 


A 


i*'(A)*'(#o  r1/2 


0(A)  - 
A  -  p 


*'(*A)*'(fy)|*l/2 


*(®a) 


ifl. 


-  *(»;/) 
-  e 


(3.4) 


i^A 

where  A  =  e  ,  fi  =  e  ,  0  <  0^  .  0^  <  2*  .  In  particular,  if  A  and  /i  are 
distinct  eigenvalues  of  H,  then  0(A)  =  0(^)  =  0,  and  therefore  v^v^  =  0. 


Proof.  By  (2.20),  (2.6)  and  (2.7), 


’  Wj 

Aj 

v\  = 

A 

v2 

I  n-s 

(A 


I  A)_1z/<5  (A) 


(3.5) 


and  therefore 


IS 


VA%  =  6iha  -  rA)'1(A  -  I/,)‘1 


=  (HA  W/O)'1  E 


Kj 


j=l  (Aj  -  AHAj  -  (i) 


Now 


AAj  _  _  _w_*l  _  a  \ 

(Aj-A)  (Aj-//)  ”  "  A-ziUj-A  Aj-zJ 


A;  A; 


(Aj-AXAj-zi) 


(3.6) 


*  (3.7) 


Substituting  (3.7)  into  (3.6)  yields 


VA  VA< 


=  (2*(A)*(/0) 


-l 


A 

A+zi 


.  n 

2  E 
v  j  =  l 


I  <i  t  2aJ 

A-  Aj 


n 

-  2  £ 
1=1 


L£j 


J 1  AJ) 

£  —  Aj  J 


and  by  (2.19),  (3.1)  and  (3.3), 


vAMV/i  =  (^'(A)^(zi)A/i)'1'^A 


1/2,  ^ ( A)  -  <KjO 


A  -  /£ 


=  (*'(0A)$'(0/1))_l/2  i e'^A  *(*a)- 


J*A 

e  -  e 


This  shows  (3.4).  □ 


The  denominator  |A-zi|  in  (3.4)  suggests  that  it  may  be  numerically 
difficult  to  obtain  orthogonal  eigenvectors  when  the  associated 
eigenvalues  are  very  close.  The  following  lemma  sheds  some  light  on  this 
situation,  and  shows  that  due  to  deflation  the  roots  of  $(0)  are  bounded 
away  from  each  other. 

iff. 

Lemma  3 . 2  Let  Aj  =  e  ,  1  <  j  <  n,  be  the  eigenvalues  of  A,  and  let  z  = 
C<j3 j*=1  be  defined  by  (2.7).  Let  be  an  arbitrary  but  fixed  positive 
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constant  and  assume  that  the  Aj  are  pairwise  distinct  and  that  2  |  Cj  I  >  (i 
for  all  j.  These  conditions  can  be  made  valid  by  deflation.  Assume  that 
the  Aj  are  sorted  so  that  0  <  <  02  <  ■  ■  •  <  6n  <  2n ,  and  let  #n+1  :=  2k  + 

0l.  Let  9  6  [0,2tt[  be  a  zero  of  $(#),  and  let  k  be  such  that  0k  <  6  <  0k  +  1 . 
Then 


9-9k  >  Jc^k  +  I-A)  =>  °k+l~9  >  min{Ti(0k  +  l-ek)  >  7T>  * 


16  ^  k  +  l  >  3 


0k  +  \~0  >  ^(0k  +  i~0k)  =>  6-$k  >  min{I^(ek  +  1-0k)  ,  ^>  . 


.(3.8) 

(3.9) 


Proof.  Introduce  the  index  sets 


II  :=  {j  :  0  <  ffj  +  2  kI  <  6  +  r ,  for  some  l  6  Z,  1  <  j  <  n}  , 


I2  :=  {  j  s  0- *  <  9 j  +  2  k t  <  9 ,  for  some  £  £  Z,  1  <  j  <  n} 
Then  Ij  n  I2  =  0  and  I j  U  I2  =  {l,2,...,n}.  Fu rt he r 


A  +  A 


A  -  A 


-j  =  cotf-EEj 


e  ~  0..  {  <  o  ,  j  e  ix, 


V  2 


>  0  , 


j  6  I2. 


In  particular,  k  G  I2  and,  provided  that  k  <  n,  k+1  e  Ij.  If  k  =  n  then 
1  €  Ij.  Moreover, 


, 0-6 


cot 


r-2^)  *  -t(vj) . 

C-V')  »  • 


cot 


From  $(0)  =  0,  it  follows  that 


Vj  €  ll 


Vj  €  I; 


(3.10) 


-  E  I  Cj  1 2  cot(-2-J)  =  E  I  Cj  1 2  c°t(-2J) 
jeii  v  ’  jei2  w  ' 


(3.11) 
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By  (3.10) ,  (3.11)  and  z  z  =  1  we  obtain,  provided  that  k  <  n, 


-Kk+i!2  cot( — ^  co't{~Ti)  ’  (3-12) 

or,  equivalently, 

|Ck  +  il2  tan (-^)  <  tan(-^ — )  .  (3.13) 

If  k  =  n  then  we  define  (n+i  :=  and  (3 . 12) - (3 . 13)  remain  valid. 

Now  assume  that 

0  -  K  >  J(Vm  -  »k)  •  (3-14) 

We  wish  to  determine  a  lower  bound  for  0k_j_j-0.  From  tan  ^  —  x  ^or 

0  <  x  <  ^  it  follows  that  if  0  <  #*  +  1“^  <  ,  then 

tan(-k-tL---)  <  0k  +  1  -  e  .  (3.15) 

Substituting  (3 . 1 4 ) - (3 . 1 5 )  into  (3.13)  yields 

Kk  +  ll2  tan  ( 4  (^k  +  l  _  ®k )  )  ^  ^k  +  l  ~  9  » 

and  from  tan  (^  (0k_j_j-  0k))  >  ^(^k  +  l-^k)  ’  we  obtain 

I  ‘‘k+l  I  2  \  (^k  +  1  -  0k)  ^  ^k  +  l  "  9  •  (3.16) 

Finally,  substituting  l<k+1|  >  ^  into  (3.16)  yields  (3.8). 

In  order  to  show  (3.9),  we  note  that  from  (3 . 10) - (3 . 1 1 )  and  zHz  =  1  it 
fol lows  that 

-  cot( - 2*  +  l)  -  I  *>k  I  2  cot(g  2 

or,  equivalently, 

tan(- 2-^)  >  |(k  | 2  tan(<k+12~  *)  ,  (3.17) 
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which  corresponds  to  (3.13).  We  now  assume  that 


*k+l  ~  *  >  J(*k+i  -  *k)  •  (3-18) 

We  would  like  to  determine  a  lower  bound  for  9-9k  .  Similarly  as  in  the 
derivation  of  (3.15)  we  obtain  that  if  0  <  9-9k  <  ^ ,  then 

tanf - ^ — ~  j  <  9  -  9k  .  (3.19) 

F  rom  (3.17)-(3.19)  and  tan(i(#k_j_j-0))  >  Fj(k+l-^)  »  we  obt  a  i  n 

e  -  9k  >  kk|2i(*k+i  -  ek)  .  (3.20) 

Finally,  substituting  I  Ck  I  -i  into  (3.20)  yields  (3.9).  □ 


Our  final  lemma  shows  that  the  computed  eigenvectors  are  close  to 
orthogonal  if  the  zeros  of  ®{9)  are  evaluated  with  sufficient  accuracy. 


Lemma  3.3  .  Let  A  =  d  i  ag  [Aa  ,  A2  ,  .  .  .  ,  An]  ,  and  let  A  ,  /j  be  computed 
approximations  of  the  distinct  roots  A,  /i  of  <j> .  Introduce  the  relative 
errors  ok  ,  0k  of  Ak-A  and  Ak-/j,  respectively,  i.e. 


Ak  -  A  -  (Ak  -  A)(l  +  ak) 
Ak  -  £  =  (Ak  “  A* )  (  1  +  /?k) 


k  =  1,2..  .  .  ,  n  . 


(3.21) 


Assume  that  for  some  constant  0  <  <  <  1,  |  e>k  |  <  f  and  |  0k  |  <  c  for  all  k, 

and  that  |  A  |  =  ( /i  (  =  1  .  Then 


|vaHCv^|  <  f  (2  +  o({-^)2 


where  C  =  d  i  ag  [pj  ,  p7  .  .  .  .  ,  />n]  with 
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Qk  ±  4  ±  °k^k  f <HA)<H/0 

(1  +  5k)(l  +  0k)[6(X)S(p) 


(3.22) 


For  r}  :=  e'0T),  0  <6^  <  2n ,  we  define  6(  i;)  :  =  (^d>,(v)T7)1/'2  =  (5  (^>7)  )l/2  • 

Proof.  We  first  note  that  since  A,  g  are  computed  by  determining  zeros  of 
$(0),  0  <  6  <  2tt  ,  the  requirement  |A|  =  | /}  |  =  1  is  satisfied.  Analogously 
to  (3.C)  we  obtain 


v  - 

A 


(«5(A)<5(A)  )_1  E 

k=l 


Kk  1 2 

(Ak  “  A)  ( Ak  "  £) 


=  (#(A  )#(/i))-1  E 
k=l 


_ kki! _ 

(Ak  -  A)  (Ak  -  /i)(l  +  5k)(l  +  0k ) 


where  the  last  equality  follows  from  (3.21). 

0  =  v  %  =  (6(X)f(tl)y1  E  ^ - 

A  k=l  (Ak  -  A)(Ak  - 


Now 


/0 


and  (2.19)  imply  that 

f  1  Ck  1 2 

k=l  (Ak  -  A)(Ak  -  /0 


0  . 


The  ref  ore 


ICkl 


V-  V-  =  (6(\)6(ii)  )'J  E  -« - = - 

A  **  ^  k  =  l  (Ak-A)(Ak-/0(l+5k)(l+/?k) 


-  E  K* 


k=i  (Ak-AKAk-/0 


=  (A(A)^  (A)  E 


Kkl 


k=i  (Ak_AKAk-/0  h  1+f5k)  (1+&k) 


((l+n.H 


-  1 


)  ■ 


(3.23) 
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L-J 

which  shows  that  v-  v- 

A  r 


with  C  defined  by  (3.22).  F rom  ( 2 . 19)  , 


Cv 


(3.3),  (3.2)  and  (3.21)  it  follows  that 


KA)\2 

6(\)J 


Aa>a 

Aa)a 


Kk  1 2 


(-Ak)A 

(A-Ak)2 


( -Ak)  A 
(A-Ak)2 


A/A  \ 

( 1+Qk) 2 


(3 


F  rom 


( _ Ak A)  /  ( A- Ak ) 2 


A/A  > 

(1  +  Qk)2  - 


(1 


>  0  it  follows  that  A/(A(l+ak)2)  >  0  and  therefore 

+  Kl)*2  >  (1  +  O'2  •  (3. 


Substituting  (3.25)  into  (3.24)  yields  6(A)/<5(A)  <  1  +  c,  and  similarly' 
can  show  that  f>  (ft) /S  (ft)  <  1  +  f.  Hence,  by  (3.22), 

!  Pk  I  <  (1  +  02  =  f  (2  +  0(H~f)2  •  (3 

Final  1 y , 

|  v-  V-  I  =  |  Cv^l  <  II V  ^  1 1 2  II C  ||  2  || 1 1 2  =  II  ^  II 2  =  I  I  ’ 

and  the  desired  bound  now  follows  from  (3.26).  0 


.  24) 


25) 


one 


.26) 


4.  The  orthogonal  eigenproblem 

The  computational  work  required  for  the  real  orthogonal  eigenproblem 
is  smaller  than  for  the  unitary  one.  This  section  discusses  these 
differences,  and  considers  some  details  of  our  implementation  of  a  DC 
scheme  for  the  real  orthogonal  eigenproblem.  Our  computer  program  is  for 
the  case  when  H  6  Rnxn  with  n  =  2^,  where  £  is  a  positive  integer,  and  we 
assume  in  this  section  that  n  is  of  this  form.  Many  of  our  comments  are 
valid  for  more  general  values  of  n,  also. 

We  first  note  that  the  subdivision  of  the  eigenproblem  for  H  into 
smaller  e i gen prob 1 ems ,  as  described  by  (1.3)-(1.7),  does  not  require  any 
computational  work.  Subdivision  yields  the  block-diagonal  matrix 

H  :=  G^Gs.  .  .G^G^G^Gn  ,  (4.1) 


and  we  obtain  simple  formulas  for  the  eigenpairs  of  each  2x2  block  on  the 
diagonal  as  follows.  Let 


G  :  = 


-7 

<T 


<7 

7 


€  R2x2 


-1  <7<  1,  (T  >  0  ,  72  +  <t2=  1 


(4.2) 


Since  G  is  real,  symmetric ,  orthogonal  and 
{Ai,A2},  we  have  Ax  =  1  and  A2  =  -1.  Let  Xj 
unit  length.  Then  we  can  choose 

r£i  =  *2-l/2o  +  7yl/2 
U2  =  2  l/2(l  +  7 ) 1/2 , 

and  from  a  -  ( 1  -  7 2 )1^2  it  follows  that 

r£i  =  2-l/2(i  -  7) 1/2 
U2  =  «t2-i/2(1  -  7)-l/2  . 


has  distinct  eigenvalues 
=  [^x'^2]T  be  an  eigenvector  of 


(4.3) 


(4.4) 
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Cancellation  of  significant  digits  is  avoided  by  using  (4.3)  if  7  >  0  and 
(4.4)  otherwise.  An  eigenvector  associated  with  A2  =  -1  is  given  by  x2  := 


If  7 n  =  -1  then  Gn  =  In,  and  the  eigenpairs  of  G^Gn  are  those  of  Gn_! 
If  7^=1  we  need  to  determine  the  eigenpairs  of 


G' 


T'n-l  ^n-l 
<Tn-l  ~7n-l 


We  find  that  the  eigenvalue  =  -7n_i  +  i <rn. t  of  G1  has  an  associated 
eigenvector  Xj  :=  [2-1^2  ,  -2'1^2]  T  ,  and  the  eigenvalue  A2  =  -7^  -  i^n-i  *las  an 

associated  eigenvector  x2  :=  [2-1^2 , 2’1/,2j  T  . 

Note  that  since  the  eigenvalues  of  G  given  by  (4.2)  are  A  =  ±1, 
independent  of  -1  <  7  <  1,  deflation  takes  place  numerous  times  during  the 
com put at  ions. 


We  turn  to  the  computation  of  the  Householder  transformation  (1.4). 

In  oder  to  avoid  cancellation  of  significant  digits,  we  compute  {ws,u>s+1} 
given  by  (1.6)  as  follows.  If  7S  >  0,  then  we  use  (1.6a)  and  replace 
( 1 . 6b)  by 

ws  +  i  :=  -  crs21^2(l  +  7S)  1 ^2  .  (1.6b*) 

In  case  7S  <  0,  we  use  (1.6b)  and  replace  (1.6a)  by 

ws  :=  <ts2'i/2(1  -  7s)'l/2  •  (1.6a') 


Due  to  H  having  real-valued 
of  H  occur  in  complex  conjugate 
6  <  8  <  7r  have  to  be  computed  . 

*(*)  =  £  Kj  I  cot(^-2—) 


elements,  the  eigenvalues  and  eigenvectors 
pairs.  Therefore  only  zeros  of  $(0)  for 
Mo  reove  r , 
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can  be  simplified.  Assume  that  0  <  0k  <  it  for  some  k  <  n  and  let 


^k  +  l  -  #|< 


Then  |  <k  |  =  |  Ck -f-i  I  >  anc*  we  obtain 


"cot 


(V)  +  l<k  +  ,l2cot(^ 


:)  =  k»i’(<»*(V)  -  c»t(5f)) 


+8  \ 


_  i  /*  |  2  2  sin  9  _  |  _ 

'  cos  6  -  cos  6 1.  *<  .  ?0\,+9\  .  /6i ,-0\ 


sin 


sin(— )sinl~!r) 


(4.5) 


We  use  the  right  hand  side  of  (4.5)  in  the  computations. 

we  need  to  evaluate  cot^^)  as  well  as  cot|^j  =  tan 

The  contribution  from  (4.5)  to  $'(&)  is 

91,  ,2  d(  sin  6  \  =  o\(  '9  1  ~  cos  6  COS  ^ 

'  d^lcos  8  -  cos  8V)  'k 


(cos  8  -  cos  6 „)2 


If  0k  =  0  then 


(4.6) 


The  stable  evaluation  of  the  right  hand  side  of  (4.6)  can  be  accomplished 
as  described  in  Table  4.1. 


Cond i t i ons _ 

cck  <  0 

cck  >  0  and  g-^g;  >  0 
cck  >  0  and  s+s:  <  0 


Evaluate 


Tab  1 e  4.1: 
c  :=  cos  0, 


Stable  evaluation 
ck  :  =  cos  0k ,  s  :  = 


of  (l-cck)/(c-ck)2 , 
sin  9 ,  sk  :=  sin  8 k 


where 


The  interlacing  of  the  zeros  of  < £(A)  with  the  impiies  that  it 

easily  can  be  determined  whether  8=0  or  8  =  n  are  zeros  of  $(0) •  Let 
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"the  6^  ,  1  <  k  <  n ,  be  ordered  so  that  0  <  61  <  $2  <  •  •  •  <  dp  <  *  < 


7P+l 


<  ...  <  6n  <  .  Since  the  Ak  =  exp(i#k)  appear  in  complex  conjugate 

pairs,  we  obtain 


(6X  > 

0  => 

*(0) 

-  0  , 

Up  < 

v  => 

*(») 

=  0  . 

Final ly ,  we 

consider  the 

computation  of 

by (2 . 20) . 

Let 

=: 

cw<‘». 

w(1)l 

w^ 

J 

e  cs  , 

W2  =: 

[wf, 

(2), 

■  »  wn-s  J  » 

w<2> 

J 

e  cn‘s  , 

s„(l) 


:.(!)- 


„(1) 


Ax  =:  d  iag  [exp(  )  ,exp(  )  ,  .  .  .  ,exp(i^g  )]  ,  0  <  9^  <  2x  , 


a(2) 


A2  =:  diag[exp(i0j  )  ,exp(i02  )  »  •  •  •  >exp(i0n-s)]  i  0  <  0j  <  2n  , 

W1Hesws  =:  K(I1),41),...,<i1)JT  , 

y  H  r,(2)  ,(2)  (2)  T 

"2  elws+i  -•  L^i  »  •  •  •  ><.n-sJ  » 

and  A  =:  exp(i<?)  ,  0  <  6  <  2?r .  Then 


W1(I-A1HA)'1W1He 


=  e  (i  -  exp (i(^-^j1))))1 

j=l  J  J  J 


,0)  a) 


E  (1  -  exp(i(M}1))))-1({1)*j1) 

+  E  [(l-exp(i(fl-ff|1))))-1<j1)wj1)  +  (l-exp(i(fl+ej1))))-1<jiyi)] 

/  1  \  J  JJ  J  J  J 

O<0j  ;<7T 
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1  E  Cl  +  icot(|))C]  V]1'  +  i  E  Cl-i^n(|))C  ,;wp 

<f=o  %}!U  J 


(I)  (i)  sin  0- 

+  £  [(*««'  M  ’)  +  - 77T — 1 - m— 

0<4')Or  „  .  A  +<\  .  A 

J  2  s  i  n( -“-g — Isinf-^ — J 


WJ  In(<i  ci  )] 


Al)+e, 


4  si"  •  £  (si"!^— )sin(J“5-))'‘  Re«jI,“j‘))  • 


We  may  assume  that  close  eigenvalues  have  been  eliminated  from  and  A 

by  deflation,  and  that  therefore  the  and  $j2^  are  distinct.  Hence,  th 

sums  over  ^  =  0  and  ^  =  t r  contain  at  most  one  term  each. 

Analogously  to  (4.7)  we  obtain 

W2(A2-IA)'1W2He1o,s+1  =  "f  (exp(  i^|2)) -exp(  i6)  )-1c[2)w[2) 

j_l  j  j  j 


=  J  L  (l  +  icot(|))42)wJ2)  +  J  E  (-l  +  itan(|))c|2)w5 

C91  ^  J  J  ^  Z.  J  J 


<f=o 


0<»!2>CV  . 

J  S  1  1 


,+,f) 

1  s  i  n| 

W2> 

f  J 

,  2  J 

i 

i,  2 

J  1 

r 

1  s  i  n| 

f-r 

,  2  J 

l  2 

hr)  Re(cf)wj‘ 


(4.8 


ie'*  sin  ^|2)(s  1  n(  t- ) s  1  n(~ir~ )) 1 

O<0j  <JT 


+  2  s  i n 


0<e!2)<5T 


».#<2>  »  «(2) 


The  evaluation  of  <5(A)  defined  by  (2.19)  can  also  be  simplified.  We  have 


6(X) 


'<j 


(i) 


|  A -exp  (  i  )| 


n-s 

(Ti —  +  ^ 

j  =  l 


«<j 


(2) 


1 1/2 


|  A-exp(  i^j2^)  I  2' 


where ,  e . g . , 


8 

£ 


1C 


(i) 


=  l  £  1C: 1}  I  2  /  sin2  (^)  +  \  £  K]lj|2/cos2(|) 


(1) 


j=l  |A-exp(i^(1))  |2 


*fU 


(1) 


(4.9) 


+  i  e  icfV 
0«f<* 


6  a(1)  e+fl(1) 

(sin(~j_ )) 2  -  (^2M-r 


The  simplifications  of  this  section  for  the  orthogonal  eigenproblem  have 
been  implemented  in  a  Pascal  program.  Several  other  mathematically 
equivalent  forms  of  (4.7)-(4.9)  could  also  be  used.  We  have  tried  to  find 
formulas  that  avoid  unnecessary  loss  of  significant  digits. 
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5. 


Numerical  examples 


We  report,  results  of  some  computed  examples  with  an  experimental 
program  for  the  orthogonal  e i genprob 1 em .  The  program  is  written  in  Turbo 
Pascal  4.0  and  was  run  on  an  IBM  PC  AT  computer  with  unit  roundoff 
u  =  2‘39  ss  2-10'12.  Our  code  implements  the  formulas  of  Section  4. 

Generally  very  accurate  answers  are  obtained.  Lemma  3.2  indicates, 
however,  that  a  zero  9  of  $(9)  may  be  very  close  to  a  singular  point  9 j  of 
$(0)  and  by  Lemma  3.3  the  difference  9-6-  has  to  be  computed  to  high 
re  1  at i ve  accuracy  in  order  to  yield  nearly  orthogonal  eigenvectors. 

Example  5.2  below  shows  that,  indeed,  9-0-  can  be  extremely  tiny  and  that 
loss  of  accuracy  in  both  eigenvectors  and  eigenvalues  may  result.  This 
loss  of  accuracy  could  be  reduced,  e.g.  ,  by  representing  9  and  9j  in  higher 
precision  arithmetic. 

In  this  section  A  G  CnXx  denotes  the  diagonal  matrix  with  the  computed 
eigenvalues  of  H  G  RnXn  as  entries,  and  W  G  CnXn  is  the  matrix  with  the 
computed  eigenvectors.  We  evaluate  the  residual  errors  HHW'-WAHoq  and 
||  WH  W  —  1 1|  oc  ,  where  ||  Hqq  denotes  the  uniform  matrix  norm. 

Example  5.1.  This  example  discusses  the  application  of  the  unitary  and 
orthogonal  e i gen prob 1 ems  to  the  construction  of  Gauss-Szego  quadrature 
rules.  Consider  the  inner  product  on  the  unit  circle 


<f  ,g> 


I 


|A|  =  1 


f(A)  g(A)  da  (A) 


(5.1) 


with  a  positive  measure  da(A).  Let  t/-k  ,  0  <  k  <  n,  be  mon  i  c  orthogonal 
polynomials  with  respect  to  (5.1).  They  satisfy  a  recurrence  relation 
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V’o(A)  :=  1 

V’k(A)  :=  AV'k_i(A)  +  7kV’k-i(A),  1  <  k  <  n, 


(5.1a) 


(5.2b) 


for  some  parameters  7k  €  C  such  that  ( 7k  |  <  1  for  1  <  k  <  n.  Here 

V’  k_  i  ( A )  :=  Ak'1V’k.1  ( 1/A)  is  the  ’’reversed  polynomial.”  Let  7n  €  C  be  an 

arbitrary  complex  number  of  unit  magnitude,  and  define  t^n  by  (5.2b)  with 
k  :=  n.  Writing  the  recursions  (5.2)  (for  1  <  k  <  n)  in  matrix  form 
yields  the  unitary  matrix 

H  =  GiG2  •  .  . GnlGn ,  (5.3) 

whose  eigenvalues  {Ak}k_j  are  the  zeros  of  ^<n  •  Here  Gk  is  defined  by  7k 

according  to  (1.2)  for  1  <  k  <  n.  Hence,  the  parameters  {7k}n  are  the 

Schur  parmaeters  for  H.  Let  H  =  WAWH  be  a  spectral  resolution,  and  define 
the  weights  p k  :=  |e^Wek|  for  1  <  k  <  n.  Then 

f  f  ( A)  da  (A)  =  £  Pk  f(Ak)  +  «n(f) 

•*  |  A  |= 1  k=l 

is  a  Gauss-Szego  quadrature  rule  with  respect  to  the  measure  da(A), 
because  the  error  <n(T)  vanishes  when  f  is  any  trigonometric  polynomial  of 
degree  less  than  n.  See  [Gr2]  for  details.  The  computed  examples 

illustrate  the  case  when  all  Schur  parameters  7k  are  real  valued  and  H 

therefore  is  real  orthogonal. 

A  particularly  simple  example  i s  7k  :=  0,  1  <  k  <  n,  and  yn  :=  -1. 

Then  V’k(A)  =  Ak  ,  0  <  k  <  n ,  and  V’n(A)  =  An  -  1,  and  therefore 

f Ak  =  exp(2iri  (k-l)/n)  , 

<  1  <  k  <  n  .  (5.4) 

lpk  =  1/n  . 
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These  Schur  parameters  have  been  used  for  Table  5.1.  In  the  table  “# 
defl.  close  e.v.”  stands  for  number  of  deflations  due  to  close 


eigenvalues,  and  “#  defl.  small  |(k|”  is  short  for  number  of  deflations 
due  to  components  (k  of  z  of  small  magnitude.  Two  eigenvalues  are 
considered  close  if  (2.26)  is  satisfied  for  :=  1-10  ,  and  I  (k  I  i® 

regarded  small  if  (2.24)  is  valid  for  :=  1-10'5.  These  values  of  and 
are  used  in  all  computed  examples  of  this  section. 


n 

IIHW-WAHoo 

l|WHW-I||0o 

#  defl.  close  e.v. 

#  def 1 .  smal 1 

4 

4.6- 10‘12 

1 . 5  - 10’ 11 

2 

0 

8 

6.41 0"12 

3 . 0  - 10" 11 

8 

0 

16 

1 . 4- 10’ 11 

5 . 4  - 10" 11 

24 

0 

32 

2.910’11 

1  . 8- 10'10 

64 

0 

64 

3.910’11 

3.4-10’10 

160 

0 

Table  5.1:  7k  :=  0,  1  <  k  <  n;  -)rn  :=  -1 

For  7k  :=  0,  1  <  k  <  n,  and  7n  :=  1,  we  obtain  the  polynomials  V’k(A)  = 

Ak  ,  0  <  k  <  n  ,  and  ^n(A)  :  =  An  +  1.  Hence,  the  eigenvalues  are  Ak  = 

exp(  i7r(2k-l)/n)  ,  1  <  k  <  n,  and  the  Gauss-Szego  weights  pk  are  the  same  as 

in  (5.4).  Table  5.2  shows  computations  for  the  present  Schur  parameters, 
and  differs  from  Table  5.1  mainly  in  that  fewer  deflations  take  place. 


n 

IIHW-WAHoo 

f|WHW-  I  ||oo 

#  defl.  close  e.v.  #  defl. 

smal  1 

4 

7.8-10’12 

1 . 610"11 

0 

0 

8 

1 .71 0’11 

4 , 2-1 0'11 

2 

0 

16 

3.  MO’11 

1 .610’10 

10 

0 

32 

4.  M0'11 

3.51  O’10 

34 

0 

64 

5 . 6  - 10' 11 

7.5-10’10 

98 

0 

Tab 1 e  5.2: 

1“H 

II 

c 

i- 

C 

V 

VI 

T— < 

o 

u 
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In  Tables  5. 3-5. 4  we  have  chosen  7 k  :=  0.8,  1  <  k  <  n.  This  makes  the 


Ak  gathe 

r  in  the 

left  half  plane. 

For  the  examples  of  Table 

5 . 3  we 

max  Re 

Ak  <  -1- 

For  the  exampl 

es 

of  Table  5.4  we  obtain  max 

Re  Ak 

^k  ^  1 

K  4 

n 

IIHW-WAIIoo 

||WHW-I||oo 

# 

defl.  close  e.v.  #  defl. 

smal  1 

4 

2.7-1  O'12 

1 . 6-10'11 

2 

0 

8 

5.5-10'11 

1  . 8- 10‘10 

8 

0 

16 

5.2-10'11 

3 . 2- 10'10 

24 

0 

32 

3 . 2  -1  O'8 

9 . 3- 10'8 

63 

1 

64 

3.2-108 

1 . 6- 10'7 

157 

3 

Table  5.3:  7k 

:  = 

0.8,  1  <  k  <  n  ;  7n  :=-l 

n 

IIHW-WAIIoo 

I|WHW-I||00 

# 

defl.  close  e.v.  #  defl. 

smal  1 

4 

4 . 8- 10'11 

1 .7-1  O'10 

0 

0 

8 

9 . 4  - 10' 11 

5.5-10'10 

2 

0 

16 

4.8-1  O'10 

6.6-1 0'9 

10 

0 

32 

6 .3-10  10 

2.3-10'8 

34 

0 

64 

4 . 2  - 10~8 

1 . 9  -1 0'7 

97 

1 

Table  5.4:  7k 

;  - 

:  0.8,  1  <  k  <  n ;  7n  :=1 

Ki. 


In  the  last  computed  quadrature  rules  of  this  example  we  let  the  7k  ,  1 

<  k  <  n,  be  uniformly  distributed  in  the  open  interval  ]-l,l[,  and  let  7n 
be  -1  or  1  with  probability  ^  each.  The  7k  are  determined  with  the  random 
number  generator  of  Pascal.  Table  5.5  shows  the  result  of  30 
e igenprobl ems  so  generated.  The  maximum,  average  and  minimum  in  Table  5.5 
are  over  all  30  e igenproblems . 
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Wi-1 


IIHW-WAHoo 

l|WHW-I||00 

#  defl.  close  e.v. 

#  def 1  .  smal 1  |  ( k 

max 

7.210'7 

2.5-10'6 

30 

0 

ave  rage 

5 . 8  - 10'8 

2 .310“7 

26.5 

0 

min 

2.91 0‘9 

1 . 5- 10'8 

22 

0 

Table  5.5:  Uniformly  distributed  7k  €]-l,l[,  1  <  k  <  n;  uniformly 

distributed  7n  6  {-1,1}.  Max,  average  and  min  are  over 
30  e i genprob 1 ems  with  n  :=  32 

The  numerical  experiments  of  Table  5.5  indicate  that  for  many  choices 
of  Schur  parameters  7k  ,  the  magnitudes  |£k|  a-re  not  sufficiently  small  to 
give  rise  to  frequent  deflations.  This  behavior  has  also  been  observed  in 
many  other  computed  experiments.  In  contrast,  massive  deflation  in  DC 
methods  for  symmetric  tridiagonal  matrices  often  is  caused  by  small 
components  of  the  vector  correspnding  to  z  =  □ 


Example  5.2.  This  example  suggests  that  it  might  not  be  possible  to 


i nc  rease 

:  the  smal 1 

lower  bound 

for  min  |  0-0,  | 

i  J 

of  Lemma  3 . 2 

s ign if i cant ly . 

The  Schu 

ir  parameters  for  Table 

5.6  are  obtained  by  reversing  the  sign  of 

the  7k  , 

1  <  k  <  n  , 

of  Table  5.4 

• 

min  |  0k  | 

#  defl . 

#  defl  . 

n 

1  <k<n 

IIHW-WAHoo 

l|WHW-I||oo 

close  e.v. 

small  |Ck  | 

4 

6.61  O'2 

6 . 9  - 10“ 11 

3.11 0“10 

0 

0 

8 

8 . 1-10"4 

7 . 2  ■  1 0~10 

2 . 5'  1 0'8 

2 

0 

16 

0* 

1 . 2- 10‘7 

3.11 0“8 

10 

0 

32 

0* 

7.210'7 

1 .9-1 0“6 

34 

2 

64 

0* 

7.21 0"7 

2 . 6' 1 0"6 

97 

5 

Table  5.6:  7k  :=  -0.8,  1  <  k  <  n;  7n  :=  1.  *The  matrix  has 

numerically  the  eigenvalue  A  =  1  of  multiplicity  two. 
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Because  crk  =  0.6  >0,  1  <  k  <  n,  the  matrix  H  has  distinct  eigenvalues 

mathematically.  Numerically  two  eigenvalues  are  so  close  that  they  are 
not  distinguished  with  our  present  choice  of  (2  =  1-10’5.  A  smaller  value 
of  t2  ,  such  as  «2  =  1-10‘6,  gave  in  some  numerical  experiments  larger 
residual  errors  HHW-WAHqo  or  ||WHW-I||oc.  □ 
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